####function that computes estimation error given simulated data and target moments
function Estimation_error(data_simul::SharedMatrix{Float64}, moments::Vector{Matrix{Float64}})    
    error = 0 #preallocate
    nmoments = 582
    err_vec = zeros(nmoments)
    err_ind = 1

    #########age profiles of wages, LFP, employment, broken up by educational status#############
    age_profiles_emp, age_profiles_income = moments[2], moments[3]
    
    #loop over ages we care about: 23-65
    for t = 23:65
        restrict_t = data_simul[data_simul[:,5] .== t,:] #simulated folks of given age
        restrict_t_hs = restrict_t[restrict_t[:,2] .== 1,:] #high schoolers
        restrict_t_sc = restrict_t[restrict_t[:,2] .== 2,:] #some college
        restrict_t_coll = restrict_t[restrict_t[:,2] .== 3,:] #college

        restrict_lf_hs = restrict_t_hs[restrict_t_hs[:, 14] .== 1, :] #high schoolers in LF
        restrict_lf_sc = restrict_t_sc[restrict_t_sc[:, 14] .== 1, :] #some college in LF
        restrict_lf_coll = restrict_t_coll[restrict_t_coll[:, 14] .== 1, :] #college in LF

        restrict_emp_hs = restrict_t_hs[restrict_t_hs[:,15] .== 1, :] #employed high schoolers
        restrict_emp_sc = restrict_t_sc[restrict_t_sc[:,15] .== 1, :] #employed some college
        restrict_emp_coll = restrict_t_coll[restrict_t_coll[:,15] .== 1, :] #employed college
        
        row_moment_data = t-22 #what we'll be using in comparison

        ####compute error of LFP
        nilf_mom_hs, nilf_mom_sc, nilf_mom_coll = age_profiles_emp[row_moment_data, 4], age_profiles_emp[row_moment_data, 7], age_profiles_emp[row_moment_data, 10]
        nilf_mom_hs_s = 1 - mean(restrict_t_hs[:,14])
        nilf_mom_sc_s = 1 - mean(restrict_t_sc[:,14])
        nilf_mom_coll_s = 1 - mean(restrict_t_coll[:,14])

        #update error vector
        err_vec[err_ind] = (nilf_mom_hs_s - nilf_mom_hs) / nilf_mom_hs
        err_ind += 1
        err_vec[err_ind] = (nilf_mom_sc_s - nilf_mom_sc) / nilf_mom_sc
        err_ind += 1
        err_vec[err_ind] = (nilf_mom_coll_s - nilf_mom_coll) / nilf_mom_coll
        err_ind += 1

        ####compute error of unemployment. We were goofing this up before.
        unemp_mom_hs = (age_profiles_emp[row_moment_data, 3] / (1 - age_profiles_emp[row_moment_data, 4]))
        unemp_mom_sc = (age_profiles_emp[row_moment_data, 6] / (1 - age_profiles_emp[row_moment_data, 7]))
        unemp_mom_coll = (age_profiles_emp[row_moment_data, 9] / (1 - age_profiles_emp[row_moment_data, 10]))

        unemp_mom_hs_s = 1 - mean(restrict_lf_hs[:,15])
        unemp_mom_sc_s = 1 - mean(restrict_lf_sc[:,15])
        unemp_mom_coll_s = 1 - mean(restrict_lf_coll[:,15])

        #update error vector
        err_vec[err_ind] = (unemp_mom_hs_s - unemp_mom_hs) / unemp_mom_hs
        err_ind += 1
        err_vec[err_ind] = (unemp_mom_sc_s - unemp_mom_sc) / unemp_mom_sc
        err_ind += 1
        err_vec[err_ind] = (unemp_mom_coll_s - unemp_mom_coll) / unemp_mom_coll
        err_ind += 1

        ####compute error of wages
        w_mom_hs, w_mom_sc, w_mom_coll = age_profiles_income[row_moment_data, 2], age_profiles_income[row_moment_data, 4], age_profiles_income[row_moment_data, 6]
        w_mom_hs_s = mean(restrict_emp_hs[:,16])
        w_mom_sc_s = mean(restrict_emp_sc[:,16])
        w_mom_coll_s = mean(restrict_emp_coll[:,16])

        #update error vector
        err_vec[err_ind] = (w_mom_hs_s - w_mom_hs) / w_mom_hs
        err_ind += 1
        err_vec[err_ind] = (w_mom_sc_s - w_mom_sc) / w_mom_sc
        err_ind += 1
        err_vec[err_ind] = (w_mom_coll_s - w_mom_coll) / w_mom_coll
        err_ind += 1
    end

    println("1:", sum(err_vec.^2)) #sum up!)

    ##########wage profiles over ability and college attainment, from NLSY###########
    w_ratios_nlsy = moments[10]
    restrict_emp = data_simul[data_simul[:,15].==1, :] #earnings data for employed people
    for s = 1:3 #loop over schooling

        #get NLSY wage moments
        ratio = w_ratios_nlsy[s]
        
        ###restrict simulated data to desired age range
        restrict_s = restrict_emp[restrict_emp[:, 2] .==s, :]
        restrict_t_1 = restrict_s[restrict_s[:, 5] .<= 34, :]
        restrict_t = restrict_t_1[restrict_t_1[:, 5] .>= 24, :]

        #get low and high ability slices
        restrict_θ_low = restrict_t[restrict_t[:, 1] .<0.33, :]
        restrict_θ_high = restrict_t[restrict_t[:, 1] .>0.66, :]

        #compute ratio 
        ratio_s = mean(restrict_θ_high[:,16]) / mean(restrict_θ_low[:,16])

        #update error vector
        err_vec[err_ind] = ((ratio_s - ratio) / ratio) * 5
        err_ind += 1
    end

    println("2:", sum(err_vec.^2)) #sum up!)


    ###########college attendance probabilities over ability, income and wealth############
    #define parent income/wealth quintiles
    attendance_overall = moments[5]
    attendance_θ = moments[4]
    attendance_I, attendance_H = moments[8], moments[6]
    attendance_θ_I, attendance_θ_H = moments[9], moments[7]

    I_quartiles = quantile(data_simul[:,7], collect(0.25:0.25:1.0))
    H_quartiles = quantile(data_simul[:,8], collect(0.25:0.25:1.0))
    nθ, ns = 3, 5

    #focus on good data entries when necessary
    data_simul_22 = data_simul[data_simul[:, 5].==22, :]

    #attendance overall
    for s = 1:ns 
        restrict_s = data_simul_22[data_simul_22[:,10].==s,:]
        share_s_s = length(restrict_s[:,1]) / length(data_simul_22[:, 1])
        share_s_nlsy = attendance_overall[s]
        err_vec[err_ind] = ((share_s_s - share_s_nlsy)/share_s_nlsy)*2
        err_ind += 1
    end

    #attendance over ability
    for i = 1:nθ
        θ_low = (0.33) * (i-1)
        θ_high = 0.33 * i

        restrict_θ_1 = data_simul_22[data_simul_22[:, 1] .>= θ_low, :]
        restrict_θ = restrict_θ_1[restrict_θ_1[:, 1] .<= θ_high, :]

        for s = 1:ns
            restrict_s = restrict_θ[restrict_θ[:,10].==s,:]
            share_s_s = length(restrict_s[:,1]) / length(restrict_θ[:, 1])
            share_s_nsly = attendance_θ[i,s+1] 

            nan_check = 1 * (share_s_nsly>0) + 1*(share_s_nsly<0) 
            if nan_check == 0
                share_s_nsly = 0
            end

            err_vec[err_ind] = (share_s_s - share_s_nsly)/share_s_nsly
            err_ind += 1
        end
    end

    println("3:", sum(err_vec.^2)) #sum up!)


    #attendance over quartiles
    
    for j = 1:4
        restrict_I = data_simul_22[data_simul_22[:,7].<=I_quartiles[1],:] #construct data with parents in j'th income/house val quint
        restrict_H = data_simul_22[data_simul_22[:,8].<=H_quartiles[1],:]

        if j>1
            restrict_I_pt1 = data_simul_22[data_simul_22[:,7].<=I_quartiles[j],:]
            restrict_H_pt1 = data_simul_22[data_simul_22[:,8].<=H_quartiles[j],:]
            restrict_I = restrict_I_pt1[restrict_I_pt1[:,7].>I_quartiles[j-1],:]
            restrict_H = restrict_H_pt1[restrict_H_pt1[:,8].>H_quartiles[j-1],:]
        end

        for s = 1:ns
            #income
            restrict_s_I = restrict_I[restrict_I[:,10].==s,:]
            share_s_I = length(restrict_s_I[:,1]) / length(restrict_I[:, 1])

            nan_check = 1 * (share_s_I>0) + 1*(share_s_I<0) 
            if nan_check == 0
                share_s_I = 0
            end

            share_s_nsly_I = attendance_I[j,s+1]
            err_vec[err_ind] = (share_s_I - share_s_nsly_I) / share_s_nsly_I
            err_ind += 1
           
            #wealth
            restrict_s_H = restrict_H[restrict_H[:,10].==s,:]
            share_s_H = length(restrict_s_H[:,1]) / length(restrict_H[:, 1])

            nan_check = 1 * (share_s_H>0) + 1*(share_s_H<0) 
            if nan_check == 0
                share_s_H = 0
            end

            share_s_nsly_H = attendance_H[j,s+1]
            err_vec[err_ind] = (share_s_H - share_s_nsly_H) / share_s_nsly_H
            err_ind += 1
        end
    end
    
    println("4:", sum(err_vec.^2)) #sum up!)

    
    #attendance over ability AND quartiles jointly. A bit less weight for these, since there are a lot of them!
    for i = 1:3, j = 1:4
        θ_low = (0.33) * (i-1)
        θ_high = 0.33 * i

        restrict_θ_1 = data_simul_22[data_simul_22[:, 1] .>= θ_low, :]
        restrict_θ = restrict_θ_1[restrict_θ_1[:, 1] .<= θ_high, :]
        
        restrict_I = restrict_θ[restrict_θ[:,7].<=I_quartiles[1],:] #construct data with parents in j'th income/house val quint
        restrict_H = restrict_θ[restrict_θ[:,8].<=H_quartiles[1],:] #construct data with parents in j'th income/house val quint

        if j>1
            restrict_I_pt1 = restrict_θ[restrict_θ[:,7].<=I_quartiles[j],:]
            restrict_H_pt1 = restrict_θ[restrict_θ[:,8].<=H_quartiles[j],:]
            restrict_I = restrict_I_pt1[restrict_I_pt1[:,7].>I_quartiles[j-1],:]
            restrict_H = restrict_H_pt1[restrict_H_pt1[:,8].>H_quartiles[j-1],:]
        end

        for s = 1:ns
            
            #income/ability
            restrict_s_I = restrict_I[restrict_I[:,10].==s,:]
            share_s_I = length(restrict_s_I[:,1]) / length(restrict_I[:, 1])

            nan_check = 1 * (share_s_I>0) + 1*(share_s_I<0) 
            if nan_check == 0
                share_s_I = 0
            end

            share_s_nsly_I = attendance_θ_I[(i-1)*4 + j,s+2]
            err_vec[err_ind] = ((share_s_I - share_s_nsly_I) / share_s_nsly_I) * (1/3)
            err_ind += 1

            #wealth/ability
            restrict_s_H = restrict_H[restrict_H[:,10].==s,:]
            share_s_H = length(restrict_s_H[:,1]) / length(restrict_H[:, 1])

            nan_check = 1 * (share_s_H>0) + 1*(share_s_H<0) 
            if nan_check == 0
                share_s_H = 0
            end

            share_s_nsly_H = attendance_θ_H[(i-1)*4 + j,s+2]
            err_vec[err_ind] = ((share_s_H - share_s_nsly_H) / share_s_nsly_H) * (1/3)
            err_ind += 1
        end
    end
    
    println("5:", sum(err_vec.^2)) #sum up!)


    ###############Additional Moments -- Give a bit more weight for these guys, since they identify specific parameters###############
    moments_extra = moments[1]
    data_simul_clean = data_simul[data_simul[:, 2].!=0, :] #clean data

    ########labor force stickiness
    #switch rate
    data_clean_23 = data_simul_clean[data_simul_clean[:, 5].>22, :]
    data_clean_60 = data_clean_23[data_clean_23[:, 5].<61, :]

    switch_rate_s = mean(data_clean_60[:, 20])
    err_vec[err_ind] = ((switch_rate_s - moments_extra[1]) / moments_extra[1]) * 5
    err_ind += 1

    #employment persistence
    restrict_prev_emp = data_simul[data_simul[:, 6].== 2, :] #employed last period
    restrict_prev_nemp =  data_simul[data_simul[:, 6].!= 2, :] #not employed last period
    
    restrict_l = restrict_prev_emp[restrict_prev_emp[:, 14] .== 1, :] #supplying labor
    restrict_l_nemp = restrict_prev_nemp[restrict_prev_nemp[:, 14] .== 1, :] #supplying labor
    emp_frac = (1 - mean(restrict_l[:, 15])) / (1 - mean(restrict_l_nemp[:, 15]))

    err_vec[err_ind] = ((emp_frac - moments_extra[2]) / moments_extra[2]) * 5
    err_ind += 1

    ########Lovenheim moments
    #1
    restrict_Lovenheim = data_simul_22[data_simul_22[:,8].>0,:]
    elast_love_1 = mean(restrict_Lovenheim[:, 13])
    err_vec[err_ind] = ((elast_love_1 - moments_extra[3]) / moments_extra[3]) * 5
    err_ind += 1

    ###2
    restrict_70 = restrict_Lovenheim[restrict_Lovenheim[:,7].<1.75, :] #1.75 = 70k
    elast_love_2 = mean(restrict_70[:, 13])
    err_vec[err_ind] = ((elast_love_2 - moments_extra[4]) / moments_extra[4]) * 2 #less weight due to bigger standard error in paper
    err_ind += 1

    #########Hilger moments
    elast_hilg = mean(data_simul_22[:, 12]) #because target is positive
    err_vec[err_ind] = ((elast_hilg - moments_extra[5]) / moments_extra[5]) * 5
    err_ind += 1

    println("6:", sum(err_vec.^2)) #sum up!)

    #########Debt Moments
    data_simul_25 = data_simul[data_simul[:, 5].==25, :]

    #mean debt
    debt_mean = mean(data_simul_25[:, 24])
    err_vec[err_ind] = ((moments_extra[6] - debt_mean)/moments_extra[6]) * 2    
    err_ind += 1

    #mean debt, conditional on having debt
    data_simul_25_n0 = data_simul_25[data_simul_25[:, 24] .> 0, :]
    debt_mean_n0 = mean(data_simul_25_n0[:, 24])
    err_vec[err_ind] = ((moments_extra[7] - debt_mean_n0)/moments_extra[7]) * 2    
    err_ind += 1

    #share with debt 
    share_debt = length(data_simul_25_n0[:, 1]) / length(data_simul_25[:, 1])
    err_vec[err_ind] = ((moments_extra[8] - share_debt)/moments_extra[8]) * 2    
    err_ind += 1

    #share with debt by parent income
    for j = 1:4
        restrict_I = data_simul_25[data_simul_25[:,7].<=I_quartiles[1],:] #construct data with parents in j'th income/house val quint
        
        if j>1
            restrict_I_pt1 = data_simul_25[data_simul_25[:,7].<=I_quartiles[j],:]
            restrict_I = restrict_I_pt1[restrict_I_pt1[:,7].>I_quartiles[j-1],:]
        end

        restrict_debt = restrict_I[restrict_I[:, 24] .> 0, :] 
        share_debt_I = length(restrict_debt[:, 1]) / length(restrict_I[:, 1])
        err_vec[err_ind] = ((moments_extra[8+j] - share_debt_I)/moments_extra[8+j]) 
        err_ind += 1
    end

    println("7:", sum(err_vec.^2)) #sum up!)

    #amass and return
    error = sum(err_vec.^2) #sum up!
    error, err_vec #return estimation error
end

#####Objective function for optimization routine
function Objective_function(guess::Array{Float64,1}, utilities::Vector{Matrix{Float64}}, moments::Vector{Matrix{Float64}}; 
    ret::Int64 = 0, rec::Int64 = 0, cfact::Int64 = 0, pbr::Int64 = 1, nsim::Int64 = 2500, dlim::Int64 = 0, decomp::Int64 = 0, rec_pol::Int64 = 0)
    error = 0

    #report how we're doing
    println("")
    println("Current Guess")
    println(round.(guess, digits = 5))

    #go for it. First some restrictions on parameter guesses to prevent weird stuff.
    if guess[1]<0 || guess[2]<0 || guess[6]<=0 || guess[24]<=0 
        error = Inf
    else
        println("Solving model . . .")
        prim, est, res = Solve_model(guess, utilities; rec = rec, cfact = cfact, pbr = pbr, dlim=dlim, decomp=decomp, rec_pol=rec_pol) #solve the models

        println("Simulating . . .")
        data_simul = Simulate(prim, est, res, utilities, nsim; rec = rec, cfact = cfact, dlim=dlim, decomp=decomp, rec_pol=rec_pol) #simulation
        
        println("Computing error . . .")
        error, err_vec = Estimation_error(data_simul, moments) #estimation error

        #garbage collection
        finalize(res.val_func_emp)
        finalize(res.val_func_unemp)
        finalize(res.val_func_nilf)
        finalize(res.pol_func_emp)
        finalize(res.pol_func_unemp)
        finalize(res.pol_func_nilf)
        finalize(res.pol_func_p_I)
        finalize(res.pol_func_p_W)
        finalize(res.pol_func_p_τ)
        @everywhere GC.gc()
    end


    println("Current error: ", round(error, digits = 5))
    println("")


    #return deliverables
    deliverables = error
    if ret == 1
        deliverables = [error, err_vec, data_simul] #return more stuff
    end
    deliverables #return
end

#########
